/*
PhyML:  a program that computes maximum likelihood phylogenies from
DNA or AA homologous sequences.

Copyright (C) Stephane Guindon. Oct 2003 onward.

All parts of the source except where indicated are distributed under
the GNU public licence. See http://www.opensource.org for details.

 */

#include <unistd.h>
#include <getopt.h>
#include "utilities.h"
#include "options.h"
#include "cl.h"
#include "models.h"
#include "free.h"
#include "interface.h"
#include "mg.h"
#include "m4.h"


/*********************************************************/
/**
* Fill the Option fields, with the argc array
*/
void Read_Command_Line(option *io, int argc, char **argv)
{
	int c;
	int open_ps_file;
	int use_gamma;
	int writemode;
	int nbl;
	nbl = 0;		 // how many BL categories did the user specify? 0 = no specs given at CLI

	struct option longopts[] =
	{
			{"n_rgrft",           required_argument,NULL,0},
			{"n_globl",           required_argument,NULL,1},
			{"max_dist",          required_argument,NULL,2},
			{"n_optim",           required_argument,NULL,3},
			{"n_best",            required_argument,NULL,4},
			{"model",             required_argument,NULL,5},
			{"search",            required_argument,NULL,6},
			{"datatype",          required_argument,NULL,7},
			{"multiple",          required_argument,NULL,8},
			{"input",             required_argument,NULL,9},
			{"support",           required_argument,NULL,10},
			{"ts/tv",             required_argument,NULL,11},
			{"nclasses",          required_argument,NULL,12},
			{"pinv",              required_argument,NULL,13},
			{"alpha",             required_argument,NULL,14},
			{"inputtree",         required_argument,NULL,15},
			{"min_diff_lk_local", required_argument,NULL,16},
			{"min_diff_lk_global",required_argument,NULL,17},
			{"steph_spr",         no_argument,NULL,18},
			{"brent_it_max",      required_argument,NULL,19},
			{"rand_start",        no_argument,NULL,20},
			{"n_rand_starts",     required_argument,NULL,21},
			{"sequential",        no_argument,NULL,22},
			{"inside_opt",        no_argument,NULL,23},
			{"p_moves",           required_argument,NULL,24},
			{"fast_nni",          no_argument,NULL,25},
			{"g_pars",            no_argument,NULL,26},
			{"r_seed",            required_argument,NULL,27},
			{"collapse_boot",     required_argument,NULL,28},
			{"random_boot",       required_argument,NULL,29},
			{"print_trace",       no_argument,NULL,30},
			{"print_site_lnl",    no_argument,NULL,31},
			{"cov",               no_argument,NULL,32},
			{"cov_delta",         required_argument,NULL,33},
			{"cov_alpha",         required_argument,NULL,34},
			{"cov_ncats",         required_argument,NULL,35},
			{"ps",                no_argument,NULL,36},
			{"cov_free",          no_argument,NULL,37},
			{"no_gap",            no_argument,NULL,38},
			{"n_rr_branch",       required_argument,NULL,39},
			{"append",            no_argument,NULL,40},
			{"no_five_branch",    no_argument,NULL,41},
			{"pars_thresh",       required_argument,NULL,42},
			{"min_diff_lk_move",  required_argument,NULL,43},
			{"hybrid",            no_argument,NULL,44},
			{"use_median",        no_argument,NULL,45},
			{"run_id",            required_argument,NULL,46},
			{"pars",              no_argument,NULL,47},
			{"quiet",             no_argument,NULL,48},
			{"measure",           no_argument,NULL,49}, // --measure will print *a lot* of explicit data about the PhyML run; think of it as being hyper-verbose
			{"blclasses",         required_argument,NULL,51},
            {"blprops",           required_argument,NULL,94},
			{"opt_blprops",       required_argument,NULL,93},
			{"max_alpha",         required_argument,NULL,58},
			{"modeltest", 		required_argument,NULL,200},
			{"modeltest_criterion", 		required_argument,NULL,201},

			{"opt_algorithm",	  	required_argument,NULL,96},
			{"opt_bl",				required_argument,NULL,1008},
			{"opt_topo",			required_argument,NULL,1009},
			{"opt_pinvar",			required_argument,NULL,1010},
			{"opt_rr",				required_argument,NULL,1011},
			{"opt_state_freq",		required_argument,NULL,1012},
			{"opt_kappa",			required_argument,NULL,1017},
			{"opt_lambda",			required_argument,NULL,1018},
			{"force_NNI_before_blopt", required_argument,NULL,1020},



			{"bfgs_itsmax",			required_argument,NULL,1013},
			{"print_opt_trace",		no_argument,NULL,1014},
			{"bfgs_stepsize",		required_argument,NULL,1015},
			{"lnsearch_stepmax",    required_argument,NULL,1016},
			{"topo_opt_full",	    required_argument,NULL,1019}, // == 1, then will optimize all free params after branch swap, rather than local-branch only optimization.

			{"sa_iters_per_stage",   required_argument,NULL,52},
			{"sa_num_anneal_stages", required_argument,NULL,53},
			{"sa_scale_temp",         required_argument,NULL,54},
			{"sa_temp_end",          required_argument,NULL,55},
			{"sa_set_back",          required_argument,NULL,56},
			{"sa_temp_start",        required_argument,NULL,57},
			{"sa_sigma_brlen",       required_argument,NULL,59},
			{"sa_sigma_pinvar",      required_argument,NULL,60},
			{"sa_sigma_gamma",       required_argument,NULL,61},
			{"sa_sigma_kappa",       required_argument,NULL,87},
			{"sa_sigma_lambda",       required_argument,NULL,88},
			{"sa_sigma_blprop",			required_argument,NULL,89},
			{"sa_sigma_pi",			required_argument,NULL,90},
			{"sa_sigma_rr",			required_argument,NULL,91},
			{"sa_sigma_emig",        required_argument,NULL,62},
			{"sa_prob_nni",          required_argument,NULL,63},
			{"sa_prob_spr",          required_argument,NULL,64},
			{"sa_prob_brlen",        required_argument,NULL,65},
			{"sa_prob_gamma",        required_argument,NULL,66},
			{"sa_prob_kappa",        required_argument,NULL,67},
			{"sa_prob_lambda",       required_argument,NULL,68},
			{"sa_prob_rr",           required_argument,NULL,69},
			{"sa_prob_blprop",		required_argument,NULL,70},
			{"sa_prob_topology",     required_argument,NULL,71},
			{"sa_prob_pinvar",       required_argument,NULL,72},
			{"sa_prob_pi",           required_argument,NULL,73},
			{"sa_prob_emig",         required_argument,NULL,74},
			{"sa_tau_start",		  required_argument,NULL,75},
			{"sa_tau_end",			  required_argument,NULL,76},
			{"sa_sele_bl",		  	  required_argument,NULL,79},
			{"sa_sele_blprop",		  required_argument,NULL,80},
			{"sa_sele_gamma",		  required_argument,NULL,81},
			{"sa_sele_pinvar",		  required_argument,NULL,82},
			{"sa_sele_kappa",		  required_argument,NULL,83},
			{"sa_sele_lambda",		  required_argument,NULL,84},
			{"sa_sele_rr",		  	  required_argument,NULL,85},
			{"sa_sele_pi",		  	  required_argument,NULL,86},
			{"sa_preopt",			  required_argument,NULL,92},
			{"sa_sele_max_params",	  required_argument,NULL,95}, // what is the maximum number of single parameters to perturb at each STA step?
			{"eb_n_gens",             required_argument,NULL,77},
			{"prev_eb_sample",	  	required_argument,NULL,78},


			{0,0,0,0}
	};

	writemode = 1;
	open_ps_file = 0;
	use_gamma = 0;
	while((c = getopt_long(argc,argv,"qi:d:m:b:n:t:f:zk:v:c:a:u:l:w:ho:s:x:g:ep:z",longopts,NULL)) != -1)
	{
		switch(c)
		{
		case 48 :
		{
			io->quiet = 1;
			break;
		}
		case 'p' : case 47 :
		{
			io->in_tree = 2;
			break;
		}
		case 46 :
		{
			io->append_run_ID = 1;
			strcpy(io->run_id_string,optarg);
			break;
		}
		case 45 :
		{
			io->mod->gamma_median = 1;
			break;
		}
		case 44 :
		{
			io->mod->s_opt->hybrid_thresh = 0;
			break;
		}
		case 43 :
		{
			io->mod->s_opt->min_diff_lk_move = atof(optarg);
			if(io->mod->s_opt->min_diff_lk_move < 0)
			{
				char choix;
				PhyML_Printf("\n. Min_diff_lk_move must be a double greater than 0.\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			break;
		}
		case 42 :
		{
			io->mod->s_opt->pars_thresh = (int)atoi(optarg);
			/* 	      if(io->mod->s_opt->pars_thresh < 0) */
			/* 		{ */
			/* 		  PhyML_Printf("\n. The parsimony threshold must be an integer greater than 0.\n"); */
			/* 		  PhyML_Printf("\n. Type any key to exit.\n"); */
			/* 		  Exit("\n"); */
			/* 		} */
			break;
		}
		case 41 :
		{
			io->mod->s_opt->opt_five_branch = 0;
			break;
		}
		case 40 :
		{
			writemode = 2;
			break;
		}
		case 39 :
		{
			char choix;
			io->mod->n_rr_branch = (int)atoi(optarg);
			if(io->mod->n_rr_branch < 1)
			{
				PhyML_Printf("\n. The number of classes must be an integer greater than 0.\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			break;
		}
		case 38 :
		{
			io->rm_ambigu = 1;
			break;
		}
		case 37 :
		{
			io->mod->s_opt->opt_cov_free_rates = 1;
			break;
		}
		case 36 :
		{
			open_ps_file = 1;
			break;
		}
		case 35 :
		{
#ifdef M4
			io->m4_model = YES;
			if(!io->mod->m4mod) io->mod->m4mod = (m4 *)M4_Make_Light((io->mod->datatype == NT)?(4):(20));
			io->mod->m4mod->n_h = (int)atoi(optarg);

			if(io->mod->m4mod->n_h < 1)
			{
				char choix;
				PhyML_Printf("\n. The number of classes must be greater than 0.\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
#endif
			break;
		}
		case 34 :
		{
#ifdef M4
			io->m4_model = YES;
			if(!io->mod->m4mod) io->mod->m4mod = (m4 *)M4_Make_Light((io->mod->datatype == NT)?(4):(20));
			io->mod->m4mod->use_cov_alpha = 1;
			io->mod->m4mod->use_cov_free  = 0;

			if(!strcmp(optarg,"e") || !strcmp(optarg,"E") ||
					!strcmp(optarg,"estimated") || !strcmp(optarg,"ESTIMATED"))
			{
				io->mod->s_opt->opt_cov_alpha = 1;
				io->mod->m4mod->alpha         = 1.0;
			}
			else
			{
				io->mod->m4mod->alpha = (double)atof(optarg);

				if(io->mod->m4mod->alpha < 1.E-5)
				{
					char choix;
					PhyML_Printf("\n. The value of alpha must be greater than 1.E-5.\n");
					PhyML_Printf("\n. Type any key to exit.\n");
					if(!scanf("%c",&choix)) Exit("\n");
					Exit("\n");
				}
			}
#endif
			break;
		}

		case 33 :
		{
#ifdef M4
			io->m4_model = YES;
			if(!io->mod->m4mod) io->mod->m4mod = (m4 *)M4_Make_Light((io->mod->datatype == NT)?(4):(20));

			if(!strcmp(optarg,"e") || !strcmp(optarg,"E") ||
					!strcmp(optarg,"estimated") || !strcmp(optarg,"ESTIMATED"))
			{
				io->mod->s_opt->opt_cov_delta = 1;
				io->mod->m4mod->delta         = 1.0;
			}
			else
			{
				io->mod->m4mod->delta = (double)atof(optarg);

				if(atof(optarg) < 1.E-10)
				{
					char choix;
					PhyML_Printf("\n. The value of delta must be larger than 1.E-10.\n");
					PhyML_Printf("\n. Type any key to exit.\n");
					if(!scanf("%c",&choix)) Exit("\n");
					Exit("\n");
				}
			}
#endif
			break;
		}
		case 32 :
		{
#ifdef M4
			io->m4_model = YES;
#endif
			break;
		}
		case 31 :
		{
			io->print_site_lnl = 1;
			break;
		}
		case 30 :
		{
			io->print_trace = 1;
			break;
		}
		case 29 :
		{
			io->random_boot_seq_order = (int)atoi(optarg);
			break;
		}
		case 28 :
		{
			io->collapse_boot = (int)atoi(optarg);
			break;
		}
		case 27 :
		{
			io->r_seed = (int)atoi(optarg);
			break;
		}
		case 26 :
		{
			io->mod->s_opt->general_pars = 1;
			break;
		}
		case 25 :
		{
			io->mod->s_opt->fast_nni = 1;
			break;
		}
		case 24 :
		{
			io->mod->s_opt->p_moves_to_examine = (double)atof(optarg);
			break;
		}
		case 23 :
		{
			io->mod->s_opt->wim_inside_opt = 1;
			break;
		}
		case 0 :
		{
			io->mod->s_opt->wim_n_rgrft = atoi(optarg);
			break;
		}
		case 1 :
		{
			io->mod->s_opt->wim_n_globl = atoi(optarg);
			break;
		}
		case 2 :
		{
			io->mod->s_opt->wim_max_dist = atoi(optarg);
			break;
		}
		case 3 :
		{
			io->mod->s_opt->wim_n_optim = atoi(optarg);
			break;
		}
		case 4 :
		{
			io->mod->s_opt->wim_n_best = atoi(optarg);
			break;
		}
		case 16 :
		{
			io->mod->s_opt->min_diff_lk_local = atof(optarg);
			break;
		}

		case 17 :
		{
			io->mod->s_opt->min_diff_lk_global = atof(optarg);
			break;
		}
		case 18 :
		{
			io->mod->s_opt->steph_spr = 0;
			io->mod->s_opt->greedy    = 1;
			break;
		}
		case 19 :
		{
			io->mod->s_opt->brent_it_max = atoi(optarg);
			break;
		}
		case 20 :
		{
			io->mod->s_opt->random_input_tree = 1;
			break;
		}
		case 21 :
		{
			io->mod->s_opt->random_input_tree = 1;
			io->mod->s_opt->n_rand_starts = atoi(optarg);
			if(io->mod->s_opt->n_rand_starts < 1) Exit("\n. Nunmber of random starting trees must be > 0.\n\n");
		}
		case 's':case 6:
		{
			if((!strcmp(optarg,"spr")) || (!strcmp(optarg,"SPR")))
			{
				io->mod->s_opt->topo_search = SPR_MOVE;
				io->user_topo 						= 1;
				io->mod->s_opt->greedy      = (io->mod->s_opt->steph_spr)?(0):(1);
			}
			else if((!strcmp(optarg,"nni")) || (!strcmp(optarg,"NNI")))
			{
				io->mod->s_opt->topo_search         = NNI_MOVE;
				io->user_topo 						= 1;
				io->mod->s_opt->random_input_tree   = 0;
			}
			else if((!strcmp(optarg,"best")) || (!strcmp(optarg,"BEST")))
			{
				io->mod->s_opt->topo_search = BEST_OF_NNI_AND_SPR;
				io->user_topo 						= 1;
				io->mod->s_opt->greedy      = (io->mod->s_opt->steph_spr)?(0):(1);
			}
			break;
		}

		case 'd':case 7:
			if(!strcmp(optarg,"nt"))
				/*         if(atoi(optarg) == NT) */
			{
				io->mod->datatype = NT;
				io->mod->stepsize = 1;
				io->mod->ns = 4;

				if(
						(io->mod->whichmodel == LG)       ||
						(io->mod->whichmodel == WAG)       ||
						(io->mod->whichmodel == DAYHOFF)   ||
						(io->mod->whichmodel == JTT)       ||
						(io->mod->whichmodel == BLOSUM62)  ||
						(io->mod->whichmodel == MTREV)     ||
						(io->mod->whichmodel == RTREV)     ||
						(io->mod->whichmodel == CPREV)     ||
						(io->mod->whichmodel == DCMUT)     ||
						(io->mod->whichmodel == VT)        ||
						(io->mod->whichmodel == MTMAM)     ||
						(io->mod->whichmodel == MTART)     ||
						(io->mod->whichmodel == HIVW)      ||
						(io->mod->whichmodel == HIVB)      ||
						(io->mod->whichmodel == CUSTOMAA)
				)
				{
					io->mod->whichmodel = HKY85;
					strcpy(io->mod->modelname, "HKY85\0");
				}
			}
			else if (!strcmp(optarg,"aa"))
				/*         else if (atoi(optarg) == AA) */
			{
				io->mod->datatype         = AA;
				io->mod->stepsize         = 1;
				io->mod->s_opt->opt_kappa = 0;
				io->mod->ns               = 20;
				if(
						(io->mod->whichmodel == JC69)   ||
						(io->mod->whichmodel == K80)    ||
						(io->mod->whichmodel == F81)    ||
						(io->mod->whichmodel == HKY85)  ||
						(io->mod->whichmodel == F84)    ||
						(io->mod->whichmodel == TN93)   ||
						(io->mod->whichmodel == GTR)    ||
						(io->mod->whichmodel == CUSTOM)
				)
				{
					io->mod->whichmodel = LG;
					strcpy(io->mod->modelname, "LG\0");
				}
			}
			else
			{
				char choix;
				PhyML_Printf("\n. Unknown argument to -d option: please use `nt' for DNA or `aa' for Amino-Acids\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}

			break;

		case 'm': case 5 :
		{
			if(!isalpha(optarg[0]))
			{
				strcpy(io->mod->custom_mod_string,optarg);

				if(strlen(io->mod->custom_mod_string) != 6)
				{
					Warn_And_Exit("\n. The string should be of length 6.\n");
				}
				else
				{
					Translate_Custom_Mod_String(io->mod);
				}

				io->mod->datatype = NT;
				io->mod->whichmodel = CUSTOM;
				strcpy(io->mod->modelname, "custom");
				io->mod->s_opt->opt_kappa     = 0;
				io->mod->s_opt->opt_rr        = 1;
				io->mod->s_opt->opt_num_param = 1;
			}

			if (strcmp(optarg, "JC69") == 0)
			{
				io->mod->datatype = NT;
				io->mod->whichmodel = JC69;
				strcpy(io->mod->modelname, "JC69");
				io->mod->s_opt->opt_kappa  = 0;
			}
			else if(strcmp(optarg, "K80") == 0)
			{
				io->mod->datatype = NT;
				io->mod->whichmodel = K80;
				strcpy(io->mod->modelname, "K80");
			}
			else if(strcmp(optarg, "F81") == 0)
			{
				io->mod->datatype = NT;
				io->mod->whichmodel = F81;
				strcpy(io->mod->modelname, "F81");
				io->mod->s_opt->opt_kappa  = 0;
			}
			else if (strcmp(optarg, "HKY85") == 0)
			{
				io->mod->datatype = NT;
				io->mod->whichmodel = HKY85;
				strcpy(io->mod->modelname, "HKY85");
			}
			else if(strcmp(optarg, "F84") == 0)
			{
				io->mod->datatype = NT;
				io->mod->whichmodel = F84;
				strcpy(io->mod->modelname, "F84");
			}
			else if (strcmp (optarg,"TN93") == 0)
			{
				io->mod->datatype = NT;
				io->mod->whichmodel = TN93;
				strcpy(io->mod->modelname, "TN93");
				if(io->mod->s_opt->opt_kappa)
				{
					io->mod->s_opt->opt_lambda = 1;
				}
			}
			else if(strcmp (optarg, "GTR") == 0)
			{
				io->mod->datatype = NT;
				io->mod->whichmodel = GTR;
				strcpy(io->mod->modelname, "GTR");
				io->mod->s_opt->opt_kappa = 0;
				io->mod->s_opt->opt_rr    = 1;
			}
			else if(strcmp(optarg, "Dayhoff") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = DAYHOFF;
				strcpy(io->mod->modelname, "Dayhoff");
			}
			else if(strcmp (optarg, "JTT") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = JTT;
				strcpy(io->mod->modelname, "JTT");
			}
			else if(strcmp(optarg, "MtREV") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = MTREV;
				strcpy(io->mod->modelname,"MtREV");
			}
			else if(strcmp (optarg, "LG") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = LG;
				strcpy(io->mod->modelname, "LG");
			}
			else if(strcmp (optarg, "WAG") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = WAG;
				strcpy(io->mod->modelname, "WAG");
			}
			else if(strcmp(optarg, "DCMut") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = DCMUT;
				strcpy(io->mod->modelname, "DCMut");
			}
			else if(strcmp (optarg, "RtREV") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = RTREV;
				strcpy(io->mod->modelname, "RtREV");
			}
			else if(strcmp(optarg, "CpREV") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = CPREV;
				strcpy(io->mod->modelname, "CpREV");
			}
			else if(strcmp(optarg, "VT") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = VT;
				strcpy(io->mod->modelname, "VT");
			}
			else if(strcmp(optarg, "Blosum62") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = BLOSUM62;
				strcpy(io->mod->modelname,"Blosum62");
			}
			else if(strcmp(optarg, "MtMam") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = MTMAM;
				strcpy(io->mod->modelname, "MtMam");
			}
			else if (strcmp(optarg,"MtArt") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = MTART;
				strcpy(io->mod->modelname, "MtArt");
			}
			else if (strcmp(optarg,"HIVw") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = HIVW;
				strcpy(io->mod->modelname, "HIVw");
			}
			else if(strcmp(optarg, "HIVb") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = HIVB;
				strcpy(io->mod->modelname, "HIVb");
			}
			else if (strcmp(optarg, "custom") == 0)
			{
				io->mod->datatype = AA;
				io->mod->whichmodel = CUSTOMAA;
				strcpy(io->mod->modelname, "Read from file");
			}
			break;
		}

		case 'a': case 14 :
		{
			use_gamma = 1;
			if ((strcmp (optarg, "e") == 0) ||
					(strcmp (optarg, "E") == 0) ||
					(strcmp (optarg, "estimated") == 0) ||
					(strcmp (optarg, "ESTIMATED") == 0))
			{
				io->mod->s_opt->opt_alpha     = 1;
				io->mod->s_opt->opt_num_param = 1;
			}
			else if ((!atof(optarg)) || (atof(optarg) < 1.E-10))
			{
				char choix;
				PhyML_Printf("\n. Alpha must be > 1.E-10.\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else
			{
				io->mod->alpha = (double)atof(optarg);
				io->mod->s_opt->opt_alpha  = 0;
			}
			break;
		}
		case 'b':case 10:
		{
			if (atoi(optarg) < -6)
			{
				char choix;
				PhyML_Printf("\n. Branch test value must be a positive integer for bootstrap, between -1 and -4 for aLRT branch test, or between -5 and -6 for posterior probability estimation.\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else
			{
				if((int)atoi(optarg) > 0)
				{
					io->ratio_test       = 0;
					io->post_probs	     = 0;
					io->mod->bootstrap   = (int)atoi(optarg);
					io->print_boot_trees = 1;

					if(io->n_data_sets > 1)
					{
						char choix;
						PhyML_Printf("\n. Bootstrap option is not allowed with multiple data sets\n");
						PhyML_Printf("\n. Type any key to exit.\n");
						if(!scanf("%c",&choix)) Exit("\n");
						Exit("\n");
					}
				}
				else if (atoi(optarg) == 0)
				{
					io->mod->bootstrap = 0;
					io->ratio_test     = 0;
					io->post_probs	   = 0;
				}
				else if(atoi(optarg) > -5)
				{
					io->mod->bootstrap = 0;
					io->ratio_test     = -(int)atoi(optarg);
					io->post_probs	   = 0;
				}
				else if (atoi(optarg) == -5)
				{
					io->post_probs	   = 2; // use an existing *.eb file
					io->mod->bootstrap = 0;
					io->ratio_test     = 0;
				}
				else if (atoi(optarg) == -6)
				{
					io->post_probs	   = 1; // use the output from the current MCMC run.
					io->mod->bootstrap = 0;
					io->ratio_test     = 0;
				}
			}
			break;
		}
		case 'c':case 12:
		{
			if ((!atoi(optarg)) || (atoi(optarg) < 0))
			{
				char choix;
				PhyML_Printf("\n. Unknown argument to -c option: the number of categories must be a positive integer\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else
			{
				io->mod->n_catg = atoi(optarg);
				if (io->mod->n_catg < 2)
				{
					io->mod->s_opt->opt_alpha = 0;
				}
			}
			break;
		}
		case 'f':
		{
			if(!strcmp(optarg,"e"))
			{
				io->mod->s_opt->opt_state_freq = 0;
			}
			else if(!strcmp(optarg,"m"))
			{
				io->mod->s_opt->opt_state_freq = 1;
			}
			else if(!isalpha(optarg[0]))
			{
				double sum;

				io->mod->s_opt->opt_state_freq  = 0;
				io->mod->s_opt->user_state_freq = 1;

				char* token;
			    /* loop until finishied */
				token = strtok(optarg, " ");
				int i = 0;
				io->mod->user_b_freq[0] = atof(token);
				while (1)
				{
						i++;
						/* extract string from string sequence */
						token = strtok(NULL, " ");

						/* check if there is nothing else to extract */
						if (token == NULL)
						{
								break;
						}

						io->mod->user_b_freq[i] = atof(token);

						/* print string after tokenized */
				}

				// next, normalize the user-specified frequencies
				// Here we operate on j residues because we don't yet know if the data is NT or AA.
				// (and thus we cannot use io->mod->ns).  Rather, we operate on residues with indices
				// j through i-1, where i was defined in the previous block as the number of
				// user-specified frequency values.
				int j = 0;
				for (j = 0; j < i; j++)
				{
					sum += io->mod->user_b_freq[j];
				}
				for (j = 0; j < i; j++)
				{
					io->mod->user_b_freq[j] /= sum;
					if (io->mod->user_b_freq[j] < .0 || io->mod->user_b_freq[j] > 1.)
					{
						Warn_And_Exit("\n. Invalid base frequencies.\n");
					}
				}
			}
			break;
		}

		case 'h':
		{
			Usage();
			break;
		}

		case 'i':case 9:
		{
			char *tmp;
			tmp = (char *) mCalloc (T_MAX_FILE, sizeof(char));
			if (strlen (optarg) > T_MAX_FILE -16)
			{
				char choix;
				strcpy (tmp, "\n. The file name'");
				strcat (tmp, optarg);
				strcat (tmp, "' is too long.\n");
				PhyML_Printf("%s",tmp);
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}

			else if (!Filexists (optarg))
			{
				char choix;
				strcpy (tmp, "\n. The file '");
				strcat (tmp, optarg);
				strcat (tmp, "' does not exist.\n");
				PhyML_Printf("%s",tmp);
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else
			{
				strcpy(io->in_seq_file, optarg);
				io->fp_in_seq = Openfile(io->in_seq_file,0);
				strcpy(io->out_tree_file,optarg);
#ifdef PHYML
				strcat(io->out_tree_file,"_phyml_tree");
#else
				strcat(io->out_tree_file,"_mc_tree.txt");
#endif
				strcpy(io->out_stats_file,optarg);
#ifdef PHYML
				strcat(io->out_stats_file,"_phyml_stats");
#else
				strcat(io->out_stats_file,"_mc_stats.txt");
#endif
			}
			Free (tmp);
			break;
		}

		case 't':case 11:
		{
			if ((io->mod->whichmodel != JC69) &&
					(io->mod->whichmodel != F81)  &&
					(io->mod->whichmodel != GTR))
			{
				if ((strcmp(optarg, "e") == 0) ||
						(strcmp(optarg, "E") == 0) ||
						(strcmp(optarg, "estimated") == 0) ||
						(strcmp(optarg, "ESTIMATED") == 0))
				{
					io->mod->kappa                 = 4.0;
					io->mod->s_opt->opt_num_param  = 1;
					io->mod->s_opt->opt_kappa      = 1;
					if (io->mod->whichmodel == TN93)
						io->mod->s_opt->opt_lambda   = 1;
				}
				else
				{
					if ((!atof(optarg)) || (atof(optarg) < .0))
					{
						char choix;
						PhyML_Printf("\n. The ts/tv ratio must be a positive number\n");
						PhyML_Printf("\n. Type any key to exit.\n");
						if(!scanf("%c",&choix)) Exit("\n");
						Exit("\n");
					}
					else
					{
						io->mod->kappa = (double)atof(optarg);
						io->mod->s_opt->opt_kappa  = 0;
						io->mod->s_opt->opt_lambda = 0;
					}
				}
			}
			break;
		}
		case 'n':case 8:
		{
			if ((!atoi(optarg)) || (atoi(optarg) < 0))
			{
				char choix;
				PhyML_Printf("\n. The number of alignments must be a positive integer\n");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else io->n_data_sets = atoi (optarg);
			break;
		}
		case 'q':case 22:
		{
			io->interleaved = 0;
			break;
		}
		case 49: // print debugging messages
		{
			io->measure = 1;
			break;
		}
		case 'l': case 51: //number of BL categories "--blclasses"
		{
			io->mod->n_l = atoi(optarg);
			if(io->mod->n_l < 1){
				PhyML_Printf("\n. The number of branch length categories (--nclasses or -l) must be at least 1.");
				PhyML_Printf("\n. Type any key to exit.\n");
				Exit("\n");
			}
			if(io->mod->n_l > MAX_BL_SET)
			{
				PhyML_Printf("\n. Your analysis is limited to the maximum number of branch length categories: %d.  You can modify this limit (#MAX_BL_SET) in the source code in utilities.h\n");
				io->mod->n_l = MAX_BL_SET;
			}
			io->mod->s_opt->opt_five_branch = 0;

			nbl = io->mod->n_l;

			break;
		}
		case 93:
		{
			int tmp = io->mod->s_opt->opt_blprops;
			io->mod->s_opt->opt_blprops = atoi(optarg);
			if(io->mod->s_opt->opt_blprops != 1 && io->mod->s_opt->opt_blprops != 0) io->mod->s_opt->opt_blprops = tmp;
			break;
		}
		case 95:
		{
			int tmp = io->sa_sele_max_params;
			io->sa_sele_max_params = atoi(optarg);
			if(io->sa_sele_max_params < 0) io->sa_sele_max_params = tmp;
			break;
		}
		case 96:
		{
			int tmp = io->opt_algorithm;
			io->opt_algorithm = atoi(optarg);
			if(io->opt_algorithm < 0) io->opt_algorithm = tmp;
			break;
		}
		case 1008:
		{
			int tmp = io->mod->s_opt->opt_bl;
			io->mod->s_opt->opt_bl = atoi(optarg);
			if (io->mod->s_opt->opt_bl < 0 || io->mod->s_opt->opt_bl > 1) io->mod->s_opt->opt_bl = tmp;
			break;
		}
		case 1009:
		{
			int tmp = io->mod->s_opt->opt_topo;
			io->mod->s_opt->opt_topo = atoi(optarg);
			if (io->mod->s_opt->opt_topo < 0 || io->mod->s_opt->opt_topo > 1) io->mod->s_opt->opt_topo = tmp;
			break;
		}
		case 1010: //pinvar
		{
			int tmp = io->mod->s_opt->opt_pinvar;
			io->mod->s_opt->opt_pinvar = atoi(optarg);
			if (io->mod->s_opt->opt_pinvar < 0 || io->mod->s_opt->opt_pinvar > 1) io->mod->s_opt->opt_pinvar = tmp;
			break;
		}
		case 1011: //rr
		{
			int tmp = io->mod->s_opt->opt_rr;
			io->mod->s_opt->opt_rr = atoi(optarg);
			if (io->mod->s_opt->opt_rr < 0 || io->mod->s_opt->opt_rr > 1) io->mod->s_opt->opt_rr = tmp;
			break;
		}
		case 1012: //pi
		{
			int tmp = io->mod->s_opt->opt_state_freq;
			io->mod->s_opt->opt_state_freq = atoi(optarg);
			if (io->mod->s_opt->opt_state_freq < 0 || io->mod->s_opt->opt_state_freq > 1) io->mod->s_opt->opt_state_freq = tmp;
			break;
		}
		case 1013: //BFGS iteration maximum
		{
			int tmp = io->mod->s_opt->bfgs_itsmax;
			io->mod->s_opt->bfgs_itsmax = atoi(optarg);
			if (io->mod->s_opt->bfgs_itsmax < 0) io->mod->s_opt->bfgs_itsmax = tmp;
			break;
		}
		case 1014: //opt_trace --- print trace information related to the optimization algorithm
		{
			io->print_opttrace = 1;
			break;
		}
		case 1015:
		{
			double tmp = io->mod->s_opt->bfgs_stepsize;
			io->mod->s_opt->bfgs_stepsize = atof(optarg);
			if (io->mod->s_opt->bfgs_stepsize < 0.0) io->mod->s_opt->bfgs_stepsize = tmp;
			break;
		}
		case 1016:
		{
			double tmp = io->mod->s_opt->lnsearch_stepmax;
			io->mod->s_opt->lnsearch_stepmax = atof(optarg);
			if (io->mod->s_opt->lnsearch_stepmax < 0.0) io->mod->s_opt->lnsearch_stepmax = tmp;
			break;
		}
		case 1017:
		{
			double tmp = io->mod->s_opt->opt_kappa;
			io->mod->s_opt->opt_kappa = atoi(optarg);
			if (io->mod->s_opt->opt_kappa < 0 || io->mod->s_opt->opt_kappa > 1) io->mod->s_opt->opt_kappa = tmp;
			break;
		}
		case 1018:
		{
			double tmp = io->mod->s_opt->opt_lambda;
			io->mod->s_opt->opt_lambda = atoi(optarg);
			if (io->mod->s_opt->opt_lambda < 0 || io->mod->s_opt->opt_lambda > 0) io->mod->s_opt->opt_lambda = tmp;
			break;
		}
		case 1019:
		{
			int tmp = io->mod->s_opt->topo_opt_full;
			io->mod->s_opt->topo_opt_full = atoi(optarg);
			if(io->mod->s_opt->topo_opt_full < 0 || io->mod->s_opt->topo_opt_full > 1) io->mod->s_opt->topo_opt_full = tmp;
			break;
		}
		case 1020:
		{
			int tmp = io->mod->s_opt->force_NNI_before_blopt;
			io->mod->s_opt->force_NNI_before_blopt = atoi(optarg);
			if(io->mod->s_opt->force_NNI_before_blopt < 0 || io->mod->s_opt->force_NNI_before_blopt > 1) io->mod->s_opt->force_NNI_before_blopt = tmp;
			break;
		}

		case 52:
		{
			int tmp = io->sa_iters_per_stage;
			io->sa_iters_per_stage = atoi(optarg);
			if(io->sa_iters_per_stage < 1) io->sa_iters_per_stage = tmp;
			break;
		}
		case 53:
		{
			int tmp = io->sa_num_anneal_stages;
			io->sa_num_anneal_stages = atoi(optarg);
			//printf("(cl.c) FOUND io->sa_num_anneal_stages = %s, %d", optarg,  io->sa_num_anneal_stages);
			if(io->sa_num_anneal_stages < 1) io->sa_num_anneal_stages = tmp;
			break;
		}
		case 54:
		{
			int tmp = io->sa_scale_temp;
			io->sa_scale_temp = atoi(optarg);
			if(io->sa_scale_temp != 0 && io->sa_scale_temp != 1) io->sa_scale_temp = tmp;
			break;
		}
		case 55:
		{
			double tmp = io->sa_temp_end;
			io->sa_temp_end = atof(optarg);
			if(io->sa_temp_end < 0.0) io->sa_temp_end = tmp; // use the default T, if the user's argument was < 0.0
			if(io->sa_temp_end == 0.0) io->sa_temp_end = 0.00000001; // the final T cannot actually be 0.0, otherwise we get error when calcuting the Boltzmann's P.
			break;
		}
		//		{"sa_temp_end",          required_argument,NULL,55},
		case 56:
		{
			int tmp = io->sa_set_back;
			io->sa_set_back = atoi(optarg);
			if(io->sa_set_back < 0) io->sa_set_back = tmp;
			break;
		}
		//		{"sa_set_back",          required_argument,NULL,56},
		case 57:
		{
			double tmp = io->sa_temp_start;
			io->sa_temp_start = atof(optarg);
			if(io->sa_temp_start < 0.0) io->sa_temp_start = tmp;
			break;
		}
		//		{"sa_temp_start",        required_argument,NULL,57},
		case 58:
		{
			double tmp = io->max_alpha;
			io->max_alpha = atof(optarg);
			if(io->max_alpha < 0.0) io->max_alpha = tmp;
			break;
		}
		//		{"max_alpha",         required_argument,NULL,58},
		case 59:
		{
			double tmp = io->sa_sigma_bl;
			io->sa_sigma_bl = atof(optarg);
			if(io->sa_sigma_bl < 0.0) io->sa_sigma_bl = tmp;
			break;
		}
		//		{"sa_sigma_bl",       required_argument,NULL,59},
		case 60:
		{
			double tmp = io->sa_sigma_pinvar;
			io->sa_sigma_pinvar = atof(optarg);
			if(io->sa_sigma_pinvar < 0.0) io->sa_sigma_pinvar = tmp;
			break;
		}
		//		{"sa_sigma_pinvar",      required_argument,NULL,60},
		case 61:
		{
			double tmp = io->sa_sigma_gamma;
			io->sa_sigma_gamma = atof(optarg);
			if(io->sa_sigma_gamma < 0.0) io->sa_sigma_gamma = tmp;
			break;

		}
		//		{"sa_sigma_gamma",       required_argument,NULL,61},
		case 62:
		{
			double tmp = io->sa_sigma_emig;
			io->sa_sigma_emig = atof(optarg);
			if(io->sa_sigma_emig < 0.0) io->sa_sigma_emig = tmp;
			break;
		}
		case 87:
		{
			double tmp = io->sa_sigma_kappa;
			io->sa_sigma_kappa = atof(optarg);
			if(io->sa_sigma_kappa < 0.0) io->sa_sigma_kappa = tmp;
			break;
		}
		case 88:
		{
			double tmp = io->sa_sigma_lambda;
			io->sa_sigma_lambda = atof(optarg);
			if(io->sa_sigma_lambda < 0.0) io->sa_sigma_lambda = tmp;
			break;
		}
		case 89:
		{
			double tmp = io->sa_sigma_blprop;
			io->sa_sigma_blprop = atof(optarg);
			if(io->sa_sigma_blprop < 0.0) io->sa_sigma_blprop = tmp;
			break;
		}
		case 90:
		{
			double tmp = io->sa_sigma_pi;
			io->sa_sigma_pi = atof(optarg);
			if(io->sa_sigma_pi < 0.0) io->sa_sigma_pi = tmp;
			break;
		}
		case 91:
		{
			double tmp = io->sa_sigma_rr;
			io->sa_sigma_rr = atof(optarg);
			if(io->sa_sigma_rr < 0.0) io->sa_sigma_rr = tmp;
			break;
		}
		case 92:
		{
			// -1 = no pre-optimization
			// 0 = pre-opt using SLM
			// 1 = pre-opt using BFGS
			int tmp = io->sa_preopt;
			io->sa_preopt = atoi(optarg);
			if(io->sa_preopt < 0 && io->sa_preopt > 2) io->sa_preopt = tmp;
			break;
		}


		//		{"sa_sigma_emig",        required_argument,NULL,62},
		case 63:
		{
			double tmp = io->sa_prob_NNI;
			io->sa_prob_NNI = atof(optarg);
			if(io->sa_prob_NNI < 0.0) io->sa_prob_NNI = tmp;
			break;
		}
		//		{"prob_nni",          required_argument,NULL,63},
		case 64:
		{
			double tmp = io->sa_prob_SPR;
			io->sa_prob_SPR = atof(optarg);
			if(io->sa_prob_SPR < 0.0) io->sa_prob_SPR = tmp;
			break;
		}
		//		{"prob_spr",          required_argument,NULL,64},
		case 65:
		{
			double tmp = io->sa_prob_brlen;
			io->sa_prob_brlen = atof(optarg);
			if(io->sa_prob_brlen < 0.0) io->sa_prob_brlen = tmp;
			break;
		}
		//		{"sa_prob_brlen",        required_argument,NULL,65},
		case 66:
		{
			double tmp = io->sa_prob_gamma;
			io->sa_prob_gamma = atof(optarg);
			if(io->sa_prob_gamma < 0.0) io->sa_prob_gamma = tmp;
			break;
		}
		//		{"sa_prob_gamma",        required_argument,NULL,66},
		case 67:
		{
			double tmp = io->sa_prob_kappa;
			io->sa_prob_kappa = atof(optarg);
			if(io->sa_prob_kappa < 0.0) io->sa_prob_kappa = tmp;
			break;

		}
		//		{"sa_prob_kappa",        required_argument,NULL,67},
		case 68:
		{
			double tmp = io->sa_prob_lambda;
			io->sa_prob_lambda = atof(optarg);
			if(io->sa_prob_lambda < 0.0) io->sa_prob_lambda = tmp;
			break;
		}
		//		{"sa_prob_lambda",       required_argument,NULL,68},
		case 69:
		{
			double tmp = io->sa_prob_rr;
			io->sa_prob_rr = atof(optarg);
			if(io->sa_prob_rr < 0.0) io->sa_prob_rr = tmp;
			break;
		}
		//		{"sa_prob_rr",           required_argument,NULL,69},
		case 70:
		{
			double tmp = io->sa_prob_blprop;
			io->sa_prob_blprop = atof(optarg);
			if(io->sa_prob_blprop < 0.0) io->sa_prob_blprop = tmp;
			break;
		}

		case 79:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_bl = x;
			}
			break;
		}
		case 80:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1 || x == 2){
				io->sa_sele_blprop = x;
			}
			break;
		}
		case 81:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_gamma = x;
			}
			break;
		}
		case 82:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_pinvar = x;
			}
			break;
		}
		case 83:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_kappa = x;
			}
			break;
		}
		case 84:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1){
				io->sa_sele_lambda = x;
			}
			break;
		}
		case 85:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1 || x == 2){
				io->sa_sele_rr = x;
			}
			break;
		}
		case 86:
		{
			int x = atoi(optarg);
			if (x == 0 || x == 1 || x == 2){
				io->sa_sele_pi = x;
			}
			break;
		}

		//		{"sa_prob_blprop",required_argument,NULL,70},
		case 71:
		{
			double tmp = io->sa_prob_topo;
			io->sa_prob_topo = atof(optarg);
			if(io->sa_prob_topo < 0.0) io->sa_prob_topo = tmp;
			break;
		}
		//		{"sa_prob_topo",     required_argument,NULL,71},
		case 72:
		{
			double tmp = io->sa_prob_pinvar;
			io->sa_prob_pinvar = atof(optarg);
			if(io->sa_prob_pinvar < 0.0) io->sa_prob_pinvar = tmp;
			break;
		}
		//		{"sa_prob_pinvar",       required_argument,NULL,72},
		case 73:
		{
			double tmp = io->sa_prob_pi;
			io->sa_prob_pi = atof(optarg);
			if(io->sa_prob_pi < 0.0) io->sa_prob_pi = tmp;
			break;
		}
		//		{"sa_prob_pi",           required_argument,NULL,73},
		case 74:
		{
			double tmp = io->prob_emig;
			io->prob_emig = atof(optarg);
			if(io->prob_emig < 0.0) io->prob_emig = tmp;
			break;
		}
		case 75:
		{
			double tmp = io->sa_tau_start;
			io->sa_tau_start = atof(optarg);
			if(io->sa_tau_start < 0.0) io->sa_tau_start = tmp;
			break;
		}
		case 76:
		{
			double tmp = io->sa_tau_end;
			io->sa_tau_end = atof(optarg);
			if(io->sa_tau_end < 0.0) io->sa_tau_end = tmp;
			break;
		}
		case 77:
		{
			int tmp = io->eb_n_gens;
			io->eb_n_gens = atoi(optarg);
			if(io->eb_n_gens < 1) io->eb_n_gens = tmp;
			break;
		}
		case 78:
		{
			char *tmp;
			tmp = (char *)mCalloc(T_MAX_FILE, sizeof(char));
			if (strlen (optarg) > T_MAX_FILE -11)
			{
				char choix;
				strcpy (tmp, "\n. The file name'");
				strcat (tmp, optarg);
				strcat (tmp, "' is too long.\n");
				PhyML_Printf("%s",tmp);
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else if (! Filexists (optarg))
			{
				char choix;
				strcpy (tmp, "\n. The file '");
				strcat (tmp, optarg);
				strcat (tmp, "' doesn't exist.\n");
				PhyML_Printf("%s",tmp);
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else
			{
				strcpy(io->in_eb_file, optarg);
				io->fp_in_eb = Openfile(io->in_eb_file,0);
			}
			Free(tmp);
			break;
		}
        case 94: //"--blprops"
        {

                // 1. Sanity check:
                if(optarg[0] != '['){
                        char choix;
                        PhyML_Printf ("\n. The proportions must be in the form [NUM1,NUM2,...,NUMN] with no spaces and surrounded by '[' and ']'.");
                        PhyML_Printf("\n. Type any key to exit.\n");
                        if(!scanf("%c",&choix)) Exit("\n");
                        Exit("\n");
                }

                // 2. determine how many BL categories are implied
                // (Although its tempting to just lookup io->mod->n_l, that parameter
                // is not guaranteed to have been parsed from the CLI yet).
                int proportions = 0; // use this as the counter for the numbr of props read in.
                int i = 0;
                for(i=1; i < strlen(optarg); i++)
                {
                        if (optarg[i] == ',')
                        {
                                proportions++;
                        }
                        if (optarg[i] == ']')
                        {
                                proportions++;
                                break;
                        }
                }

                // 2b. If the number of BL proportions != the specified number of branch length
                // categories, then something is wrong.
                // Sanity check:
                if (nbl > 0)
                {
                        if (proportions != nbl)
                        {
                                PhyML_Printf("\n. You specified %d branch length proportions using --blprops, but you also specified %d branch length categories using --blclasses.  Please ensure these values are equal!\n", proportions, nbl);
                                PhyML_Printf("\n. Type any key to exit.\n");
                                Exit("\n");
                        }
                }

                // 3. allocate memory for 'proportions' number of BL sets
                io->mod->bl_props = (double *)mCalloc(proportions, sizeof(double));

                // 4. finally, parse the BL proportions as floating-point values...
                proportions = 0; // reset counter
                i = 0;                   // reset counter
                int j = 0;
                char aprop[100];
                for(i=1; i < strlen(optarg); i++)
                {
                        if(optarg[i] != ',' && optarg[i] != ']')
                        {
                                aprop[j] = optarg[i];
                                j++;
                        }
                        else
                        {
                                aprop[j+1] = '\0'; // end the string
                                io->mod->bl_props[proportions] = atof(aprop);
                                proportions++;
                                int qqq = 0;
                                for (qqq = 0; qqq<j+1;qqq++)
                                {
                                	aprop[qqq] = '\0'; //reset aprop
                                }
                                j = 0; //reset j
                        }
                }
                // Normalize the proportions:
                double sum = 0.0;
                for(i=0; i < proportions; i++)
                {
                	if (io->mod->bl_props[i] < 0.0 ) // negative proportions are not allowed.
                	{
                		io->mod->bl_props[i] *= -1;
                	}
                	sum += io->mod->bl_props[i];
                }
                if (sum == 0.0)
                {
                    for(i=0; i < proportions; i++)
                    {
                    	io->mod->bl_props[i] = 1.0 / proportions; // set all proportions to be equal
                    }
                }
                else
                {
					for(i=0; i < proportions; i++)
					{
						io->mod->bl_props[i] = io->mod->bl_props[i] / sum;
					}
                }
                io->mod->blprops_init = 1; // the bl_props array has now been initialized.
                break;
        }
		case 'u':case 15:
		{
			char *tmp;
			tmp = (char *)mCalloc(T_MAX_FILE, sizeof(char));
			if (strlen (optarg) > T_MAX_FILE -11)
			{
				char choix;
				strcpy (tmp, "\n. The file name'");
				strcat (tmp, optarg);
				strcat (tmp, "' is too long.\n");
				PhyML_Printf("%s",tmp);
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else if (! Filexists (optarg))
			{
				char choix;
				strcpy (tmp, "\n. The file '");
				strcat (tmp, optarg);
				strcat (tmp, "' doesn't exist.\n");
				PhyML_Printf("%s",tmp);
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else
			{
				//PhyML_Printf(". debug: cl.c 1315: user tree!");
				strcpy(io->in_tree_file, optarg);
				io->in_tree = 1;
				io->fp_in_tree = Openfile(io->in_tree_file,0);
			}
			Free(tmp);
			break;
		}


		case 'v':case 13:
		{
			if ((strcmp (optarg, "e") == 0) ||
					(strcmp (optarg, "E") == 0) ||
					(strcmp (optarg, "estimated") == 0) ||
					(strcmp (optarg, "ESTIMATED") == 0)) {
				io->mod->s_opt->opt_num_param = 1;
				io->mod->s_opt->opt_pinvar    = 1;
				io->mod->invar                = 1;
			}
			else if ((atof(optarg) < 0.0) || (atof(optarg) > 1.0))
			{
				char choix;
				PhyML_Printf("\n. The proportion of invariable site must be a number between 0.0 and 1.0\n");
				PhyML_Printf("\n. Type any key to exit.");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			else
			{
				io->mod->pinvar = (double)atof(optarg);
				if (io->mod->pinvar > 0.0+MDBL_MIN)
					io->mod->invar = 1;
				else
					io->mod->invar = 0;
				io->mod->s_opt->opt_pinvar = 0;
			}
			break;
		}
		case 'o':
		{
			if(!strcmp(optarg,"tlr"))
			{
				io->mod->s_opt->opt_topo = 1;
				io->mod->s_opt->opt_bl   = 1;
			}
			else if(!strcmp(optarg,"tl"))
			{
				io->mod->s_opt->opt_topo = 1;
				io->mod->s_opt->opt_bl   = 1;
			}
			else if(!strcmp(optarg,"t"))
			{
				Warn_And_Exit("\n. You can't optimize the topology without adjusting branch length too...\n");
			}
			else if(!strcmp(optarg,"lr"))
			{
				io->mod->s_opt->opt_topo = 0;
				io->mod->s_opt->opt_bl   = 1;
			}
			else if(!strcmp(optarg,"l"))
			{
				io->mod->s_opt->opt_topo = 0;
				io->mod->s_opt->opt_bl   = 1;
			}
			else if(!strcmp(optarg,"r"))
			{
				io->mod->s_opt->opt_topo = 0;
				io->mod->s_opt->opt_bl   = 0;
			}
			else if(!strcmp(optarg,"none") || !strcmp(optarg,"n"))
			{
				io->mod->s_opt->opt_topo = 0;
				io->mod->s_opt->opt_bl   = 0;
			}
			else
			{
				char choix;
				PhyML_Printf ("\n. The optimization parameter must be 'tlr' or 'tl' or 'lr' or 'l' or 'r' or ''.");
				PhyML_Printf("\n. Type any key to exit.\n");
				if(!scanf("%c",&choix)) Exit("\n");
				Exit("\n");
			}
			break;
		}
//		case '?':
//		{
//			char choix;
//			if (isprint (optopt))
//				PhyML_Printf ("\n. Unknown option `-%c'.\n", optopt);
//			else
//				PhyML_Printf ("\n. Unknown option character `\\x%x'.\n", optopt);
//			PhyML_Printf("\n. Type any key to exit.\n");
//			if(!scanf("%c",&choix)) Exit("\n");
//			Exit("\n");
//			break;
//		}

		// modeltest
		case 200:
		{
			io->modeltest = 1;

			if(!strcmp(optarg,"hlrt") )
			{
				io->modeltest_criterion = 1;
			}
			if(!strcmp(optarg,"aic") )
			{
				io->modeltest_criterion = 0;
			}
		}

		case 201:
		{
			int tmp = atoi(optarg);
			if (tmp < 0 || tmp > 1)
			{
				PhyML_Printf("\n. Invalid valid for --modeltest_criterion.  You entered %d, but the value must be 0 or 1.  I will use AIC instead.\n", tmp);
				tmp = 0;
			}
			io->modeltest_criterion = tmp;
		}

		default:
			Usage();
		}
	}


#ifndef PHYML
	if((open_ps_file) || (io->m4_model == YES))
	{
		strcpy(io->out_ps_file,io->in_seq_file);
		strcat(io->out_ps_file, "_mc_tree.ps");
		io->fp_out_ps = Openfile(io->out_ps_file,1);
	}
#endif


	/*
	 * Here are cases and situations that must be dealt with after we've parsed
	 * the command line.  Most of the cases result in throwing an error and/or exiting.
	 */

	if (io->mod->n_l == 1)
	{
		io->mod->s_opt->opt_blprops = 0; // no need to optimize proportions if the mixture model is unused.
	}

	if (io->mod->whichmodel == GTR)
	{
		io->mod->n_diff_rr = 6;
	}

	if (io->mod->whichmodel != GTR && io->mod->datatype == NT)
	{
		io->mod->s_opt->opt_rr = 0;
	}

	if((io->mod->s_opt->n_rand_starts)           &&
			(io->mod->s_opt->topo_search == NNI_MOVE) &&
			(io->mod->s_opt->random_input_tree))
	{
		Warn_And_Exit("\n. The random starting tree option is only compatible with SPR based search options.\n");
	}

	if ((io->mod->datatype == NT) && (io->mod->whichmodel > 10))
	{
		char choix;
		PhyML_Printf("\n. Err: model incompatible with the data type. Please use JC69, K80, F81, HKY, F84, TN93 or GTR\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Warn_And_Exit("\n");
	}
	else if ((io->mod->datatype == AA) && (io->mod->whichmodel < 11))
	{
		char choix;
		PhyML_Printf("\n. Err: model incompatible with the data type. Please use LG, Dayhoff, JTT, MtREV, WAG, DCMut, RtREV, CpREV, VT, Blosum62, MtMam, MtArt, HIVw or HIVb.\n");
		PhyML_Printf("\n. Type any key to exit.\n");
		if(!scanf("%c",&choix)) Exit("\n");
		Exit("\n");
	}

	if(io->m4_model == YES)
	{
#ifdef M4
		io->mod->ns *= io->mod->m4mod->n_h;
		io->mod->use_m4mod = 1;
		M4_Make_Complete(io->mod->m4mod->n_h,
				io->mod->m4mod->n_o,
				io->mod->m4mod);
#endif
	}
	else
	{
		io->mod->s_opt->opt_cov_delta      = 0;
		io->mod->s_opt->opt_cov_alpha      = 0;
		io->mod->s_opt->opt_cov_free_rates = 0;
	}

	if (io->opt_algorithm == 2) // using STA
	{
		io->user_topo 						= 1;
	}
	if(io->opt_algorithm == 3 ) // using Empirical Bayes
	{
		io->user_topo                       = 1;
		io->mod->bootstrap 					= 0; //disable bootstrap
		io->ratio_test 						= 0; //disable aLR
		io->mod->bootstrap 					= 0; //disable bootstrap
		io->ratio_test 						= 0; //disable aLRT

		if ((io->post_probs != 1) && (io->post_probs != 2))
		{
			io->post_probs						= 1;
		}
	}
	if ( (io->post_probs == 2) && (io->fp_in_eb == NULL) )
	{
		PhyML_Printf("\n. Err. You specified to estimate posterior probabilities using an existing *.eb file, but you failed to specify which file.\n");
		PhyML_Printf("    In other words, if you specify '--support -5', you must also use the argument '--prev_eb_sample'. \n");
		PhyML_Printf("\n. Type any key to exit.\n");
		Exit("\n");
	}

	if((io->mod->s_opt->opt_cov_free_rates) && (io->mod->s_opt->opt_cov_alpha))
	{
		io->mod->s_opt->opt_cov_free_rates = 0;
		io->mod->m4mod->use_cov_alpha      = 0;
		io->mod->m4mod->use_cov_free       = 1;
	}

	if(io->print_site_lnl)
	{
		strcpy(io->out_lk_file,io->in_seq_file);
		strcat(io->out_lk_file, "_phyml_lk");
		if(io->append_run_ID) { strcat(io->out_lk_file,"_"); strcat(io->out_lk_file,io->run_id_string); }
		strcat(io->out_lk_file, ".txt");
		io->fp_out_lk = Openfile(io->out_lk_file,1);
	}

	if(io->print_trace)
	{
		strcpy(io->out_trace_file,io->in_seq_file);
		strcat(io->out_trace_file,"_phyml_trace");
		if(io->append_run_ID) { strcat(io->out_trace_file,"_"); strcat(io->out_trace_file,io->run_id_string); }
		strcat(io->out_trace_file,".txt");
		io->fp_out_trace = Openfile(io->out_trace_file,1);
	}
	if(io->print_opttrace)
	{
		strcpy(io->out_opttrace_file,io->in_seq_file);
		strcat(io->out_opttrace_file,"_phyml_opttrace");
		if(io->append_run_ID) { strcat(io->out_opttrace_file,"_"); strcat(io->out_opttrace_file,io->run_id_string); }
		strcat(io->out_opttrace_file,".txt");
		io->fp_out_opttrace = Openfile(io->out_opttrace_file,1);
	}

	if(io->opt_algorithm == 2) // STA
	{
		strcpy(io->out_annealing_file,io->in_seq_file);
		strcat(io->out_annealing_file,".annealing_trace");
		if(io->append_run_ID) { strcat(io->out_annealing_file,"_"); strcat(io->out_annealing_file,io->run_id_string); }
		strcat(io->out_annealing_file,".txt");
		io->fp_out_annealing_file = Openfile(io->out_annealing_file,1);
	}

	if(io->mod->s_opt->random_input_tree)
	{
		strcpy(io->out_trees_file,io->in_seq_file);
		strcat(io->out_trees_file,"_phyml_trees");
		if(io->append_run_ID) { strcat(io->out_trees_file,"_"); strcat(io->out_trees_file,io->run_id_string); }
		strcat(io->out_trees_file,".txt");
		io->fp_out_trees = Openfile(io->out_trees_file,1);
	}

	if((io->print_boot_trees) && (io->mod->bootstrap > 0))
	{
		strcpy(io->out_boot_tree_file,io->in_seq_file);
		strcat(io->out_boot_tree_file,"_phyml_boot_trees");
		if(io->append_run_ID) { strcat(io->out_boot_tree_file,"_"); strcat(io->out_boot_tree_file,io->run_id_string); }
		strcat(io->out_boot_tree_file,".txt");
		io->fp_out_boot_tree = Openfile(io->out_boot_tree_file,1);

		strcpy(io->out_boot_stats_file,io->in_seq_file);
		strcat(io->out_boot_stats_file,"_phyml_boot_stats");
		if(io->append_run_ID) { strcat(io->out_boot_stats_file,"_"); strcat(io->out_boot_stats_file,io->run_id_string); }
		strcat(io->out_boot_stats_file,".txt");
		io->fp_out_boot_stats = Openfile(io->out_boot_stats_file,1);
	}

	if(io->append_run_ID)
	{
		strcat(io->out_tree_file,"_");
		strcat(io->out_stats_file,"_");
		strcat(io->out_tree_file,io->run_id_string);
		strcat(io->out_stats_file,io->run_id_string);
	}
	strcat(io->out_tree_file,".txt");
	strcat(io->out_stats_file,".txt");

	// If the user enabled modeltesting using hLRT, but their data is AA,
	// then we default back to AIC instead of hLRT.
	if (io->modeltest == 1 && io->modeltest_criterion == 1 && io->mod->datatype == AA)
	{
		PhyML_Printf("\n. Warning: The hLRT modeltest is not valid for amino acid data.  I will use AIC instead of hLRT.\n");
		io->modeltest_criterion = 0;
	}

	io->fp_out_tree  = Openfile(io->out_tree_file,writemode);
	io->fp_out_stats = Openfile(io->out_stats_file,writemode);

	return;
}

/*********************************************************/
